////////////////////////////////////////////////////////////////////////
//HEADERS
////////////////////////////////////////////////////////////////////////
#include<stdio.h>
#include<stdlib.h>  
#include<string.h>  
#include<malloc.h>
#include<math.h>
#include<time.h>

//CABECERAS DE GSL
#include<gsl/gsl_math.h>
#include<gsl/gsl_errno.h>
#include<gsl/gsl_sort.h>
#include<gsl/gsl_sort_vector.h>
#include<gsl/gsl_matrix.h>
#include<gsl/gsl_spline.h>
//Minimizacion
#include <gsl/gsl_min.h>
#include <gsl/gsl_multimin.h>
#include<gsl/gsl_bspline.h>
#include <gsl/gsl_multifit.h>
#include <gsl/gsl_statistics.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_deriv.h>


////////////////////////////////////////////////////////////////////////
//MACRO
////////////////////////////////////////////////////////////////////////
#define POW gsl_pow_int
#define VERBOSE 0

#define vprintf if(VERBOSE) printf

////////////////////////////////////////////////////////////////////////
//UNITS
////////////////////////////////////////////////////////////////////////
#define CLIGHT 3E10 //cm/s
#define YEAR 365.0*24*60*60 //s
#define MSUN 1.989E33 //g
#define KPC 3.2616E3*CLIGHT*YEAR //cm
#define GYR 1E9*YEAR //s
#define KPS 1E5 //cm/s

#define UM 1E10*MSUN 
#define UL 1*KPC 
#define UV 1*KPS
#define UT 0.9785*GYR
#define GGAD 43007.1 //UL^3/(UM UT^2)
#define GCONST GGAD*UL*UL*UL/(UM*UT*UT)

////////////////////////////////////////////////////////////////////////
//DATATYPES
////////////////////////////////////////////////////////////////////////
typedef struct 
{
  int id;
  int type;
  float pos[3];
  float vel[3];
  float mass;
  float pot;
  float posc[3];
  float velc[3];
  float mang[3];
  float mr;
  float r;
  float vc;
  float epsilon;
  float Ebinding;
  float je[3];
} particulas;

typedef struct
{
  int Npart[6];
  double mass[6];
  double time;
  double   redshift;
  int      flag_sfr;
  int      flag_feedback;
  int      npartTotal[6];
  int      flag_cooling;
  int      num_files;
  double   BoxSize;
  double   Omega0;
  double   OmegaLambda;
  double   HubbleParam; 
  char     fill[256- 6*4- 6*8- 2*8- 2*4- 6*4- 2*4 - 4*8];  /* fills to 256 Bytes */
} io_header;

typedef struct
{
  int id;
  float cm[3];
  float vcm[3];
  float lcm[3];
  float M;
} CM;

typedef struct
{
  float comp[6];
}vector;

typedef struct
{
  int stellarage,metallicity,outputpotencial;
  int outputacceleration,outputchangeofentropy,outputtimestep;
}config;

typedef struct
{
  float Mtot;
  float L[3];
  float rCM[3];
  float vCM[3];
  float cos[3];
  float Ro;
}disk;

////////////////////////////////////////////////////////////////////////
//GLOBAL VARIABLES
////////////////////////////////////////////////////////////////////////
particulas **part;
particulas *particles;
io_header Header;
CM cmdummy;
vector *Pos;
vector *Vel;
int *Id;
float *Mass;
float *Softening,SofteningFactor[6];
int N_tot,N_part[6],N_min,N_max;
config Config;
disk Disk;
//double *dfdx,*d2fdx2;

size_t *indexradio,*indexaltura;  
double *radio,*altura;

int NPotential=0;
float *RPotential;
float *Potential;
float *EpicFreq;

int n0,n1,n2,n3,n4,n5;

////////////////////////////////////////////////////////////////////////
//ROUTINES
////////////////////////////////////////////////////////////////////////

//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
//INPUT / OUTPUT
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
/*
  ======================================================================
  Open a file
  ======================================================================
*/
FILE *fileOpen(char filename[],char mode[])
{
  FILE *f;

  if((f=fopen(filename,mode))==NULL){
    fprintf(stderr,"Error openning '%s' for %s\n",filename,mode);
    exit(1);
  }

  return f;
}

/*
  ======================================================================
  Create plain Header file and populate Pos,Vel,Id,masa
  
  Parameters:

  confname: name of the configuration file

  ======================================================================
*/

// organizo en orden creciente de fvector, indexando en ivector
int gsl_int_int_sort(int dimension,int *fvector,int *ivector)
{
 
  int *aux_ivector,i;
  int *aux_fvector;
  size_t *P;
 
  P=(size_t *) malloc((size_t) dimension*sizeof(size_t));
  if(P == NULL){
    printf("Allocation error routines:81\n");
    exit(0);
  }
 
  aux_fvector=(int *) malloc((size_t) dimension*sizeof(int));
  if(aux_fvector == NULL){
    printf("Allocation error routines:81\n");
    exit(0);
  }
 
  aux_ivector=(int *) malloc((size_t) dimension*sizeof(int));
  if(aux_ivector == NULL){
    printf("Allocation error routines:81\n");
    exit(0);
  }
 
  for(i=0; i<dimension; i++)
    {
      aux_ivector[i]=ivector[i];
      aux_fvector[i]=fvector[i];
    }
 
  gsl_sort_int_index(P,fvector,1,dimension);
 
  for(i=0; i<dimension; i++)
    {
      fvector[i]=aux_fvector[P[i]];
      ivector[i]=aux_ivector[P[i]];
    }
 
  free(P);
  free(aux_fvector);
  free(aux_ivector);
 
  return 0;
 
}

int readGadgetBinary(char confname[],char filename[])
{
  int i,j;
  FILE *fHeader,*fconf,*fdata;
  char tmp[1000];
  int dummy;

  ////////////////////////////////////////////////////////////////////////
  //READ CONFIGURATION FILE
  ////////////////////////////////////////////////////////////////////////
  
  fconf=fileOpen("options_makefile.dat","r");

  int stellarage,metallicity,outputpotencial;
  int outputacceleration,outputchangeofentropy,outputtimestep;

  fscanf(fconf,"%d",&stellarage);
  fscanf(fconf,"%d",&metallicity);
  fscanf(fconf,"%d",&outputpotencial);
  fscanf(fconf,"%d",&outputacceleration);
  fscanf(fconf,"%d",&outputchangeofentropy);
  fscanf(fconf,"%d",&outputtimestep);
  
  fclose(fconf);
  
  ////////////////////////////////////////////////////////////////////////
  //READ DATA FILE
  ////////////////////////////////////////////////////////////////////////
  fdata = fopen(filename,"r");
  if (fdata==NULL) printf("No se pudo abrir %s\n",filename);

  ////////////////////////////////////////////////////////////////////////
  //SAVE HEADER
  ////////////////////////////////////////////////////////////////////////
  fread(&dummy,sizeof(dummy),1,fdata);
  fread(&Header,sizeof(io_header),1,fdata);
  fread(&dummy,sizeof(dummy),1,fdata);

  sprintf(tmp,"%s_Header.dat",filename);
  printf("\n");
  printf("Reading header %s\n",tmp);
  printf("writing header in %s\n",tmp);
  printf("\n");
  fHeader=fopen(tmp,"w");
  fprintf(fHeader,"Snapshot %s\n",filename);
  fprintf(fHeader,"%d %d %d %d %d %d\n",
	  Header.Npart[0],Header.Npart[1],Header.Npart[2],
	  Header.Npart[3],Header.Npart[4],Header.Npart[5]);
  fprintf(fHeader,"Mass= %f %f %f %f %f %f\n",
	  Header.mass[0],Header.mass[1],Header.mass[2],
	  Header.mass[3],Header.mass[4],Header.mass[5]);
  fprintf(fHeader,"time %lf, redshift %lf, flag_sfr %d, flag_feedback %d\n",
	  Header.time,Header.redshift,Header.flag_sfr,Header.flag_feedback);
  fprintf(fHeader,"Npartotal= %d %d %d %d %d %d\n",
	  Header.npartTotal[0],Header.npartTotal[1],Header.npartTotal[2],
	  Header.npartTotal[3],Header.npartTotal[4],Header.npartTotal[5]);
  fprintf(fHeader,"flag_cooling %d, num_files %d, BoxSize %lf, Omega0 %lf, OmegaLambda %lf, HubbleParam %lf\n",
	  Header.flag_cooling,Header.num_files,Header.BoxSize,Header.Omega0,Header.OmegaLambda,Header.HubbleParam);
  fprintf(fHeader,"tamaño restante %lu\n",sizeof(Header.fill));
  fprintf(fHeader,"\n\n");
  fclose(fHeader);

  ////////////////////////////////////////////////////////////////////////
  //TOTAL NUMBER OF PARTICLES
  ////////////////////////////////////////////////////////////////////////
  N_tot=0;
  printf("Reading snapshot with:\n");
  for(i=0;i<6;i++)
    {
      N_part[i]=Header.npartTotal[i];
      N_tot+=N_part[i];
      printf("%d particles of type %d\n",N_part[i],i);
    }
  printf("%d particles in the snapshot\n",N_tot);

  ////////////////////////////////////////////////////////////////////////
  //ALLOCATE AND READ
  ////////////////////////////////////////////////////////////////////////
  
  particles = (particulas *)malloc((size_t)N_tot*sizeof(particulas));
  if(particles == NULL){
    printf("No se pudo alocar particles\n");
    exit(0);
  }

  //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  //POSITIONS
  //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  fread(&dummy,sizeof(dummy),1,fdata);
  for(i=0; i<N_tot; i++)
    {
      fread(&particles[i].pos[0],sizeof(float),3,fdata);
    }
  fread(&dummy,sizeof(dummy),1,fdata);

  //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  //VELOCITIES
  //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  fread(&dummy,sizeof(dummy),1,fdata);
  for(i=0; i<N_tot; i++)
    {
      fread(&particles[i].vel[0],sizeof(float),3,fdata);
    }
  fread(&dummy,sizeof(dummy),1,fdata);

  //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  //IDS
  //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  fread(&dummy,sizeof(dummy),1,fdata);
  for(i=0; i<N_tot; i++)
    {
      fread(&particles[i].id,sizeof(int),1,fdata);
    }
  fread(&dummy,sizeof(dummy),1,fdata);

  
  //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  //MASSES
  //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  fread(&dummy,sizeof(dummy),1,fdata);

  N_min=N_max=0;

  for(j=0;j<=5;j++)
    {
      N_max=N_max+Header.npartTotal[j];
      if((Header.mass[j]==0) && (Header.npartTotal[j]!=0))
	{
	  for(i=N_min;i<N_max;i++)
	    {
	      fread(&particles[i].mass,sizeof(float),1,fdata);
	    }
	}

      if((Header.mass[j]!=0) && (Header.npartTotal[j]!=0))
	{
	  for(i=N_min;i<N_max;i++)
	    {
	      particles[i].mass = Header.mass[j];
	    }
	}
      N_min=N_max;
    }
  fread(&dummy,sizeof(dummy),1,fdata);

  /*
  //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  //Potential
  //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  fread(&dummy,sizeof(dummy),1,fdata);
  for(i=0;i<N_tot;i++)
    {
      fread(&particles[i].pot,sizeof(float),1,fdata);
    }
  fread(&dummy,sizeof(dummy),1,fdata);

  fclose(fdata);
*/

  float *U,*rho,*Ne,*Nh,*h,*sfr,*Z,*stellar_age,*acce,*ecr,*timestep;
  int N_metal =  Header.npartTotal[0]+Header.npartTotal[2]+Header.npartTotal[3]+Header.npartTotal[4]+Header.npartTotal[5];
  
  U = (float *)malloc(Header.npartTotal[0]*sizeof(float));
  rho = (float *)malloc(Header.npartTotal[0]*sizeof(float));
  Ne = (float *)malloc(Header.npartTotal[0]*sizeof(float));
  Nh = (float *)malloc(Header.npartTotal[0]*sizeof(float));
  h = (float *)malloc(Header.npartTotal[0]*sizeof(float));
  sfr = (float *)malloc(Header.npartTotal[0]*sizeof(float));
  Z = (float *)malloc(N_metal*sizeof(float));
  stellar_age = (float *)malloc(Header.npartTotal[0]*sizeof(float));
  //  potential = (float *)malloc(N_tot*sizeof(float));
  acce = (float *)malloc(N_tot*sizeof(float));
  ecr = (float *)malloc(Header.npartTotal[0]*sizeof(float));
  timestep = (float *)malloc(N_tot*sizeof(float));
  
 
  //*****************************************************
  //Reading additional properties
  //*****************************************************

  if(Header.npartTotal[0]!=0)
    {
      
      //Read internal energy for particles of gas
      //====================================================
      fread(&dummy,sizeof(dummy),1,fdata);
      for(i=0;i<Header.npartTotal[0];i++)
	{
	  fread(&U[i],sizeof(float),1,fdata);
	}
      fread(&dummy,sizeof(dummy),1,fdata);
      //====================================================
      
      //Read density for particles of gas
      //====================================================
      fread(&dummy,sizeof(dummy),1,fdata);
      for(i=0;i<Header.npartTotal[0];i++)
	{
	  fread(&rho[i],sizeof(float),1,fdata);
	}
      fread(&dummy,sizeof(dummy),1,fdata);
      //====================================================
      
      if(Header.flag_cooling==1)
	{
	  //Read electron density for particles of gas
	  //====================================================
	  fread(&dummy,sizeof(dummy),1,fdata);
	  for(i=0;i<Header.npartTotal[0];i++)
	    {
	      fread(&Ne[i],sizeof(float),1,fdata);
	    }
	  fread(&dummy,sizeof(dummy),1,fdata);
	  //====================================================
	  
	  //Read neutral hydrigen density for particles of gas
	  //====================================================
	  fread(&dummy,sizeof(dummy),1,fdata);
	  for(i=0;i<Header.npartTotal[0];i++)
	    {
	      fread(&Nh[i],sizeof(float),1,fdata);
	    }
	  fread(&dummy,sizeof(dummy),1,fdata);
	  //====================================================
	}
      
      //Read SPH smoothing length for particles of gas
      //====================================================
      fread(&dummy,sizeof(dummy),1,fdata);
      for(i=0;i<Header.npartTotal[0];i++)
	{
	  fread(&h[i],sizeof(float),1,fdata);
	}
      fread(&dummy,sizeof(dummy),1,fdata);
      //====================================================
      
      if(Header.flag_sfr==1)
	{
	  //Read star formation rate for particles of gas
	  //====================================================
	  fread(&dummy,sizeof(dummy),1,fdata);
	  for(i=0;i<Header.npartTotal[0];i++)
	    {
	      fread(&sfr[i],sizeof(float),1,fdata);
	    }
	  fread(&dummy,sizeof(dummy),1,fdata);
	  //====================================================
	  
	  //Read formation time of star for particles of gas
	  //====================================================
	  if(stellarage==1)
	    {
	      fread(&dummy,sizeof(dummy),1,fdata);
	      for(i=0;i<Header.npartTotal[0];i++)
		{
		  fread(&stellar_age[i],sizeof(float),1,fdata);
		}
	      fread(&dummy,sizeof(dummy),1,fdata);
	    }
	  //===================================================
	}
      //====================================================
      
    }
  
  //Read metallicity for gas and stars
  //====================================================
  if((Header.flag_sfr==1) && (metallicity==1))
    {
      fread(&dummy,sizeof(dummy),1,fdata);
      for(i=0;i<N_metal;i++)
	{
	  fread(&Z[i],sizeof(float),1,fdata);
	}
      fread(&dummy,sizeof(dummy),1,fdata);
    }
  //====================================================
  
  //Read gravitational potential for all particles
  //====================================================
  if(outputpotencial==1)
    {
      fread(&dummy,sizeof(dummy),1,fdata);
      for(i=0;i<N_tot;i++)
	{
	  fread(&particles[i].pot,sizeof(float),1,fdata);
	}
      fread(&dummy,sizeof(dummy),1,fdata);
    }
  //====================================================
  
  //Read acceleration for all particles
  //====================================================
  if(outputacceleration==1)
    { 
      fread(&dummy,sizeof(dummy),1,fdata);
      for(i=0;i<N_tot;i++)
	{
	  fread(&acce[i],sizeof(float),1,fdata);
	}
      fread(&dummy,sizeof(dummy),1,fdata);
    }
  //====================================================
  
  //Read rate of change of entropic function for gas 
  //====================================================
  if(outputchangeofentropy==1)
    { 
      fread(&dummy,sizeof(dummy),1,fdata);
      for(i=0;i<Header.npartTotal[0];i++)
	{
	  fread(&ecr[i],sizeof(float),1,fdata);
	}
      fread(&dummy,sizeof(dummy),1,fdata);
    }
  //=====================================================
  
  //Read timestep for all particles
  //====================================================
  if(outputtimestep==1)
    { 
      fread(&dummy,sizeof(dummy),1,fdata);
      for(i=0;i<N_tot;i++)
	{
	  fread(&timestep[i],sizeof(float),1,fdata);
	}
      fread(&dummy,sizeof(dummy),1,fdata);
    }
  //====================================================
 

  //######################################
  //Sorting particles by id
  //######################################
  
  particulas *aux;
  int *ID,*Index;
  size_t N=N_tot;
  aux = (particulas *)malloc((size_t)N_tot*sizeof(particulas));
  if(aux == NULL){
    printf("Allocation error routines:81\n");
    exit(0);
  }
  Index = (int *)malloc((size_t)N_tot*sizeof(int));
  if(Index == NULL){
    printf("Allocation error routines:81\n");
    exit(0);
  }
  ID = (int *)malloc((size_t)N_tot*sizeof(int));
  if(ID == NULL){
    printf("Allocation error routines:81\n");
    exit(0);
  }

  for(i=0; i<N_tot; i++)
    {
      for(j=0; j<3; j++)
	{
	  aux[i].pos[j] = particles[i].pos[j];
	  aux[i].vel[j] = particles[i].vel[j];
	}
      aux[i].mass = particles[i].mass;
      aux[i].pot = particles[i].pot;
      aux[i].id = particles[i].id;
      ID[i] = particles[i].id;
      Index[i] = i;
    }

  gsl_int_int_sort(N,ID,Index);


  n0 = Header.npartTotal[0];
  n1 = Header.npartTotal[1];
  n2 = Header.npartTotal[2];
  n3 = Header.npartTotal[3];
  n4 = Header.npartTotal[4];
  n5 = Header.npartTotal[5];
  for(i=0; i<N_tot; i++)
    {
      for(j=0; j<3; j++)
	{
	  particles[i].pos[j] = aux[Index[i]].pos[j];
	  particles[i].vel[j] = aux[Index[i]].vel[j];
	}
      particles[i].mass = aux[Index[i]].mass;
      particles[i].pot = aux[Index[i]].pot;
      particles[i].id = aux[Index[i]].id;
      particles[i].id = i;
      
      if(i < n0) 
	particles[i].type = 0;
      
      if( (i >= n0) && ( i < (n0+n1)) ) 
	particles[i].type = 1;
      
      if( (i >= (n0+n1)) && (i < (n0+n1+n2)) ) 
	particles[i].type = 2;
      
      if( (i >= (n0+n1+n2) ) && (i < (n0+n1+n2+n3)) ) 
	particles[i].type = 3;
      
      if( (i >= (n0+n1+n2+n3)) && (i < (n0+n1+n2+n3+n4)) ) 
	particles[i].type = 4;
      
      if((i >= (n0+n1+n2+n3+n4)) && (i < (n0+n1+n2+n3+n4+n5)) ) 
	particles[i].type = 5;
    }
  
  free(Index);
  free(ID);
  free(aux);
  
  /*
    for(i=0;i<=5;i++)
    {
    printf("%d %e %e %e %e %e %e\n",
    particles[i].id,
    particles[i].pos[0],particles[i].pos[1],particles[i].pos[2],
    particles[i].vel[0],particles[i].vel[1],particles[i].vel[2]);
    }
  */
  
  printf("Done.\n\n");
  
  return 0;
}



int centerMass(int *Index, int nindex, int limlower, int limuper,CM *centerm)
{
  int i, count;
  double mtot=0.0, cm[3], vcm[3];

  for(i=0; i<3; i++)
    {
      cm[i]=0.0;
      vcm[i]=0.0;
      (*centerm).cm[i] = 0.0;
      (*centerm).vcm[i] = 0.0;
    }
  
  if(Index==NULL)
    {
      
      printf("Computing CM with interval (%d,%d)\n",limlower,limuper);
      
      mtot = 0.0;
      for(i=limlower; i<limuper; i++)
	{
	  cm[0] = cm[0] + particles[i].pos[0]*particles[i].mass;
	  cm[1] = cm[1] + particles[i].pos[1]*particles[i].mass;
	  cm[2] = cm[2] + particles[i].pos[2]*particles[i].mass;
	  
	  vcm[0] = vcm[0] + particles[i].vel[0]*particles[i].mass;
	  vcm[1] = vcm[1] + particles[i].vel[1]*particles[i].mass;
	  vcm[2] = vcm[2] + particles[i].vel[2]*particles[i].mass;
	  
	  mtot = mtot + particles[i].mass;
	}
      
      (*centerm).cm[0] = cm[0]/mtot;
      (*centerm).cm[1] = cm[1]/mtot;
      (*centerm).cm[2] = cm[2]/mtot;
      
      (*centerm).vcm[0] = vcm[0]/mtot;
      (*centerm).vcm[1] = vcm[1]/mtot;
      (*centerm).vcm[2] = vcm[2]/mtot;
      
    }
  else
    {
      
      printf("Computing CM with index\n");
      mtot = 0.0;
      count = 0;
      for(i=0; i<nindex; i++)
	{
	  cm[0] = cm[0] + particles[Index[i]].pos[0]*particles[Index[i]].mass;
	  cm[1] = cm[1] + particles[Index[i]].pos[1]*particles[Index[i]].mass;
	  cm[2] = cm[2] + particles[Index[i]].pos[2]*particles[Index[i]].mass;
	  
	  vcm[0] = vcm[0] + particles[Index[i]].vel[0]*particles[Index[i]].mass;
	  vcm[1] = vcm[1] + particles[Index[i]].vel[1]*particles[Index[i]].mass;
	  vcm[2] = vcm[2] + particles[Index[i]].vel[2]*particles[Index[i]].mass;
	  
	  mtot = mtot + particles[Index[i]].mass;
	  count++;
	}
      
      (*centerm).cm[0] = cm[0]/mtot;
      (*centerm).cm[1] = cm[1]/mtot;
      (*centerm).cm[2] = cm[2]/mtot;
      
      (*centerm).vcm[0] = vcm[0]/mtot;
      (*centerm).vcm[1] = vcm[1]/mtot;
      (*centerm).vcm[2] = vcm[2]/mtot;
      
      printf("with %d particles\n",count);
      
    }
  
  
  /*  
  printf("RCM : %e %e %e\n",(*centerm).cm[0],(*centerm).cm[1],(*centerm).cm[2]);
  printf("VCM : %e %e %e\n",(*centerm).vcm[0],(*centerm).vcm[1],(*centerm).vcm[2]);
  printf("LCM : %e %e %e\n",(*centerm).lcm[0],(*centerm).lcm[1],(*centerm).lcm[2]);
  printf("M TOTAL : %e\n\n",mtot);
  */
  
  return 0;
}

int compute_angmom(int *Index, int nindex, int limlower, int limuper,CM *centerm)
{
  
  int i,count;
  double lcm[3];
  
  for(i=0; i<3; i++)
    lcm[i]=0.0;
  
  if(Index==NULL)
    {
      
      printf("Computing LCM with interval (%d,%d)\n",limlower,limuper);
      
      
      for(i=limlower; i<limuper; i++)
	{
	  
	  lcm[0] = lcm[0] + (particles[i].pos[1]*particles[i].vel[2] 
			     - particles[i].pos[2]*particles[i].vel[1])*particles[i].mass;
	  lcm[1] = lcm[1] - (particles[i].pos[0]*particles[i].vel[2] 
			     - particles[i].pos[2]*particles[i].vel[0])*particles[i].mass;
	  lcm[2] = lcm[2] + (particles[i].pos[0]*particles[i].vel[1] 
			     - particles[i].pos[1]*particles[i].vel[0])*particles[i].mass;
	  
	}
      
      (*centerm).lcm[0] = lcm[0];
      (*centerm).lcm[1] = lcm[1];
      (*centerm).lcm[2] = lcm[2];
      
    }
  else
    {
      
      printf("Computing LCM with index\n");
      count = 0;
      for(i=0; i<nindex; i++)
	{
	  lcm[0] = lcm[0] + (particles[Index[i]].pos[1]*particles[Index[i]].vel[2] 
			     - particles[Index[i]].pos[2]*particles[Index[i]].vel[1])*particles[Index[i]].mass;
	  lcm[1] = lcm[1] - (particles[Index[i]].pos[0]*particles[Index[i]].vel[2] 
			     - particles[Index[i]].pos[2]*particles[Index[i]].vel[0])*particles[Index[i]].mass;
	  lcm[2] = lcm[2] + (particles[Index[i]].pos[0]*particles[Index[i]].vel[1] 
			     - particles[Index[i]].pos[1]*particles[Index[i]].vel[0])*particles[Index[i]].mass;
	  
	  count++;
	}
      
      (*centerm).lcm[0] = lcm[0];
      (*centerm).lcm[1] = lcm[1];
      (*centerm).lcm[2] = lcm[2];
      
      
      printf("with %d particles\n",count);
            
    }

  
  return 0;
}


int traslate(CM *zero, int *Index, int nindex)
{
  int i;

  for(i=0; i<nindex; i++)
    {
      
      particles[Index[i]].pos[0] = particles[Index[i]].pos[0] - (*zero).cm[0];
      particles[Index[i]].pos[1] = particles[Index[i]].pos[1] - (*zero).cm[1]; 
      particles[Index[i]].pos[2] = particles[Index[i]].pos[2] - (*zero).cm[2]; 
      
      particles[Index[i]].vel[0] = particles[Index[i]].vel[0] - (*zero).vcm[0];  
      particles[Index[i]].vel[1] = particles[Index[i]].vel[1] - (*zero).vcm[1];  
      particles[Index[i]].vel[2] = particles[Index[i]].vel[2] - (*zero).vcm[2];  
      
    }

  return 0;

}

//int rotate(CM *vector, int *Index, int nindex)
int rotate(float *lcm, int *Index, int nindex)
{
  
  
  int m;
  double thetar=0.0, betar=0.0, r, theta, beta;
  double M1[3][3], M2[3][3];
  float vec1[3], vec2[3], lcmr[3];
  
  
  /* en el plano y,z  --> rotacion al rededor de x*/
  
  r = sqrt( lcm[1]*lcm[1] + lcm[2]*lcm[2] );
  theta = acos(fabs(lcm[1])/r);
  theta = theta*180.0/M_PI;
    
  if((lcm[1] >= 0) && (lcm[2] >= 0))
    thetar = theta + 270.0;
  
  if((lcm[1] >= 0) && (lcm[2] < 0))
    thetar = 180.0 + (90.0 - theta);
  
  if((lcm[1] < 0) && (lcm[2] < 0))
    thetar = theta + 90.0;
  
  if((lcm[1] < 0) && (lcm[2] >= 0))
    thetar = 90.0 - theta;
  
  printf("theta=%lf thetar=%lf\n",theta,thetar);
  thetar  = thetar*M_PI/180.0;
  
  //rotate about x 
  M1[0][0] = 1.0;         M1[0][1] = 0.0;               M1[0][2] = 0.0;
  M1[1][0] = 0.0;         M1[1][1] = cos(thetar);       M1[1][2] = sin(thetar);
  M1[2][0] = 0.0;         M1[2][1] = -sin(thetar);      M1[2][2] = cos(thetar);
  
  
  for(m=0; m<nindex; m++)
    {
      
      vec1[0] = particles[Index[m]].pos[0];
      vec1[1] = particles[Index[m]].pos[1];
      vec1[2] = particles[Index[m]].pos[2];      
      
      vec2[0] = particles[Index[m]].vel[0];
      vec2[1] = particles[Index[m]].vel[1];
      vec2[2] = particles[Index[m]].vel[2];
      
      particles[Index[m]].pos[0] = M1[0][0]*vec1[0] + M1[0][1]*vec1[1] + M1[0][2]*vec1[2];
      particles[Index[m]].pos[1] = M1[1][0]*vec1[0] + M1[1][1]*vec1[1] + M1[1][2]*vec1[2];
      particles[Index[m]].pos[2] = M1[2][0]*vec1[0] + M1[2][1]*vec1[1] + M1[2][2]*vec1[2];

      particles[Index[m]].vel[0] = M1[0][0]*vec2[0] + M1[0][1]*vec2[1] + M1[0][2]*vec2[2];
      particles[Index[m]].vel[1] = M1[1][0]*vec2[0] + M1[1][1]*vec2[1] + M1[1][2]*vec2[2];
      particles[Index[m]].vel[2] = M1[2][0]*vec2[0] + M1[2][1]*vec2[1] + M1[2][2]*vec2[2]; 
      
    } 
  


  //////////////////////////////////////
  /*          Segunda rotacion        */
  //////////////////////////////////////   
  
  
  lcmr[0] = M1[0][0]*lcm[0] + M1[0][1]*lcm[1] + M1[0][2]*lcm[2];
  lcmr[1] = M1[1][0]*lcm[0] + M1[1][1]*lcm[1] + M1[1][2]*lcm[2];
  lcmr[2] = M1[2][0]*lcm[0] + M1[2][1]*lcm[1] + M1[2][2]*lcm[2];

  /* en el plano x,z --> rota al rededor de y */ 
  r = sqrt( lcmr[0]*lcmr[0] + lcmr[2]*lcmr[2] );
  beta = acos(fabs(lcmr[0])/r); 
  beta = beta*180.0/M_PI;
  
  if((lcmr[0] >= 0.0) && (lcmr[2] >= 0.0))
    betar = 90 - beta;
  
  if((lcmr[0] < 0.0) && (lcmr[2] >= 0.0))
    betar = beta + 270.0;
  
  if((lcmr[0] < 0.0) && (lcmr[2] < 0.0))
    betar = (90.0 - beta) + 180.0;
  
  if((lcmr[0] >= 0.0) && (lcmr[2] < 0.0))
    betar = beta + 90.0;
    
  printf("beta=%lf betar=%lf\n",beta,betar);
  betar = betar*M_PI/180.0;
  
  //rotate about Y
  M2[0][0] = cos(betar);      M2[0][1] = 0.0;            M2[0][2] = -sin(betar);
  M2[1][0] = 0.0;             M2[1][1] = 1.0;            M2[1][2] = 0.0;
  M2[2][0] = sin(betar);      M2[2][1] = 0.0;            M2[2][2] = cos(betar);
  

  for(m=0; m<nindex; m++)
    {
      
      vec1[0] = particles[Index[m]].pos[0];
      vec1[1] = particles[Index[m]].pos[1];
      vec1[2] = particles[Index[m]].pos[2];
      
      vec2[0] = particles[Index[m]].vel[0];
      vec2[1] = particles[Index[m]].vel[1];
      vec2[2] = particles[Index[m]].vel[2]; 
      
      particles[Index[m]].pos[0] = M2[0][0]*vec1[0] + M2[0][1]*vec1[1] + M2[0][2]*vec1[2];
      particles[Index[m]].pos[1] = M2[1][0]*vec1[0] + M2[1][1]*vec1[1] + M2[1][2]*vec1[2];
      particles[Index[m]].pos[2] = M2[2][0]*vec1[0] + M2[2][1]*vec1[1] + M2[2][2]*vec1[2];

      particles[Index[m]].vel[0] = M2[0][0]*vec2[0] + M2[0][1]*vec2[1] + M2[0][2]*vec2[2];
      particles[Index[m]].vel[1] = M2[1][0]*vec2[0] + M2[1][1]*vec2[1] + M2[1][2]*vec2[2];
      particles[Index[m]].vel[2] = M2[2][0]*vec2[0] + M2[2][1]*vec2[1] + M2[2][2]*vec2[2];
      
    } 
  

  return 0;
  
}

int particlesProperties(int numpart,int print,int time)
{
  int i,j,k,b,nbines,n,indice;
  size_t *p,*ind;
  double *r,mr,jc,v2,*E,e=0.0,Emin,Emax,J,Jmean,sigma;
  double jmax[3]={0.0,0.0,0.0},deltaE,*jz;
  double rad2grad,grad2rad;
  char namefiles[1000];

  rad2grad = 180.0/M_PI;
  grad2rad = M_PI/180.0;

  r = (double *)malloc((size_t)numpart*sizeof(double));
  p = (size_t *)malloc((size_t)numpart*sizeof(size_t));
  
  //positions
  k = b = 0;
  for(i=0; i<numpart; i++)
    {
      //R
      particles[i].posc[0] = sqrt(pow(particles[i].pos[0],2) + pow(particles[i].pos[1],2));

      //phi
      particles[i].posc[1] = atan2(particles[i].pos[1],particles[i].pos[0])*rad2grad;

      if(particles[i].posc[1]>=0.0)
	{
	  particles[i].posc[1] = particles[i].posc[1]*grad2rad;
	}
      else
	{
	  particles[i].posc[1] = (particles[i].posc[1] + 360.0)*grad2rad;
	}

      //z
      particles[i].posc[2] = particles[i].pos[2];

      //velocities
      //VR
      particles[i].velc[0] = (particles[i].pos[0]*particles[i].vel[0]+particles[i].pos[1]*particles[i].vel[1])/particles[i].posc[0];

      //Vphi
      particles[i].velc[1] = (particles[i].vel[1]*particles[i].pos[0]-particles[i].pos[1]*particles[i].vel[0])/particles[i].posc[0];

      //Vz
      particles[i].velc[2] = particles[i].vel[2];
     
      particles[i].mang[0] = (particles[i].pos[1]*particles[i].vel[2] 
			 - particles[i].pos[2]*particles[i].vel[1])*particles[i].mass;
      particles[i].mang[1] = - (particles[i].pos[0]*particles[i].vel[2] 
			 - particles[i].pos[2]*particles[i].vel[0])*particles[i].mass;
      particles[i].mang[2] = (particles[i].pos[0]*particles[i].vel[1] 
			 - particles[i].pos[1]*particles[i].vel[0])*particles[i].mass;

      r[i] = sqrt(particles[i].pos[0]*particles[i].pos[0]
		  +particles[i].pos[1]*particles[i].pos[1]
		  +particles[i].pos[2]*particles[i].pos[2]);
      particles[i].r = r[i];

      p[i] = i;
      particles[i].mr = 0.0;

      v2 = particles[i].vel[0]*particles[i].vel[0]
	+particles[i].vel[1]*particles[i].vel[1]
	+particles[i].vel[2]*particles[i].vel[2];

      particles[i].Ebinding = 0.5*particles[i].mass*v2 + particles[i].mass*particles[i].pot;
      if(particles[i].Ebinding>=0.0)
	{	
	  // printf("particulas escapando %d %lf %lf %lf\n",i,particles[i].Ebinding,0.5*particles[i].mass*v2,particles[i].mass*particles[i].pot);
	  k++;
	  if(particles[i].type!=1)
	    {
	      b++;
	    }/*
	  if(particles[i].pot>=0.0)
	    {
	       printf("OJO POTENCIAL POSITIVO\n");
	    }*/
	} 

      particles[i].epsilon = 0.0;
     
    }
  if(k>0) printf("There are %d particles with Ebinding>0.0, %d are barions\n",k,b);

  FILE *test;
  if(print==1)
    {
      sprintf(namefiles,"jz_binding_%d.dat",time);
     
      test = fopen(namefiles,"w");
    }

  E = (double *)malloc((size_t)(numpart-Header.Npart[1]-b)*sizeof(double));
  jz = (double *)malloc((size_t)(numpart-Header.Npart[1]-b)*sizeof(double));
  ind = (size_t *)malloc((size_t)(numpart-Header.Npart[1]-b)*sizeof(size_t));

  //========================
  //         NUEVO
  //========================


  //Computing Jz for types 0,2,3,4 and 5 
  j = 0;
  for(i=0; i<numpart; i++)
    {
      //if(particles[i].type!=1)
      if((particles[i].type!=1) && (particles[i].Ebinding<=0.0))
	{
	  E[j] = particles[i].Ebinding;
	  ind[j] = i;
	  jz[j] = particles[i].mang[2];
	  J = sqrt(particles[i].mang[0]*particles[i].mang[0]+
		   particles[i].mang[1]*particles[i].mang[1]+
		   particles[i].mang[2]*particles[i].mang[2]);
	  if(particles[i].type==2)
	    {
	      if(print==1)
		{
		  fprintf(test,"%e %e %e %e %e %e %e %e %lu %d\n",particles[i].pos[0],particles[i].pos[1],particles[i].pos[2],
			  particles[i].Ebinding,particles[i].mang[2],J,particles[ind[j]].Ebinding,E[j],ind[j],i);
		}
	    }
	  j++;
	}
    }


  //  printf("***particulas disco  %d : %d ***",j,numpart-Header.Npart[1]);

  //Computing m(r) and Vcir for particles
  gsl_sort_index(p,r,1,numpart);
  particles[p[0]].mr = 0.0;
  particles[p[0]].vc = 0.0;
  //particles[p[0]].epsilon = 0.0;

  /* first compute Vc */
  
  mr = particles[p[0]].mass;
  for(i=1; i<numpart; i++)
    {
      particles[p[i]].mr = mr;
      mr = mr + particles[p[i]].mass;
    }
 
  for(i=1; i<numpart; i++)
    {
      particles[p[i]].vc = sqrt(GGAD*particles[p[i]].mr/r[p[i]]);
    }
  
  for(i=1; i<numpart; i++)
    {
      jc = particles[p[i]].mass*particles[p[i]].vc*r[p[i]];
      particles[p[i]].je[2] = jc;
      //particles[p[i]].epsilon = particles[p[i]].mang[2]/jc;
    }


  //Computing Jmax=jc and epsilon for types 0,2,3,4 and 5 

  gsl_sort_index(ind,E,1,(numpart-n1-b));

  deltaE = 0.001;
  nbines = (E[ind[numpart-Header.npartTotal[1]-b-1]]-E[ind[0]])/deltaE;
  printf(" ***nbines : %d %lf %lf***\n",nbines,E[ind[0]],E[ind[numpart-Header.npartTotal[1]-b-1]]);


  FILE *test1;
  if(print==1)
    {
      sprintf(namefiles,"jcircE_%d.dat",time);
     
      test1 = fopen(namefiles,"w");
    }

  Emin = E[ind[0]];
  //printf("%d -> %lf : %lf \n",ind[0],E[ind[0]],particles[ind[0]+n0+n1].Ebinding);
  for(i=0; i<nbines; i++)
    {
      Emax = Emin + deltaE;
      // for(k=0; k<2; k++)
      //	{
      jmax[2] = -10.0;
      Jmean = 0.0;
      n = 0; 
      sigma = 0.0;
      for(j=0; j<numpart; j++)
	{
	  //if((particles[j].type!=1) && (particles[j].Ebinding<=0.0))
	  if((particles[j].type==2) && (particles[j].Ebinding<=0.0))
	    {
	      if((particles[j].Ebinding>=Emin) && (particles[j].Ebinding<Emax))
		{
		  Jmean = Jmean + particles[j].mang[2];
		  n = n+1;
		}
	    }
	}
      Jmean = Jmean/n;
      
      for(j=0; j<numpart; j++)
	{
	  //if((particles[j].type!=1) && (particles[j].Ebinding<=0.0))
	  if((particles[j].type==2) && (particles[j].Ebinding<=0.0))
	    {
	      if((particles[j].Ebinding>=Emin) && (particles[j].Ebinding<Emax))
		{
		  sigma = sigma + sqrt((particles[j].mang[2]-Jmean)*(particles[j].mang[2]-Jmean));
		}
	    }
	}
      sigma = sigma/(1.0*n-1.0);
      
      for(j=0; j<numpart; j++)
	{
	  //	  if((particles[j].type!=1) && (particles[j].Ebinding<=0.0))
	  if((particles[j].type==2) && (particles[j].Ebinding<=0.0))
	    {
	      if((particles[j].Ebinding>=Emin) && (particles[j].Ebinding<Emax))
		{
		  if((particles[j].mang[2]>jmax[2]) && (particles[j].mang[2]<=(Jmean+2.5*sigma)))
		    //if((particles[j].mang[2]>jmax[2]))
		    {
		      jmax[0] = particles[j].mang[0];
		      jmax[1] = particles[j].mang[1];
		      jmax[2] = particles[j].mang[2];
		      e = particles[j].Ebinding;
		      indice = j;
		    }
		  else
		    {
		      jmax[2] = jmax[2];
		      //  particles[j].type=8;
		    }
		}
	    }
	}
      //	}

      for(j=0; j<numpart; j++)
   	{
	  if((particles[j].Ebinding>=Emin) && (particles[j].Ebinding<Emax))
	    {
	      particles[j].epsilon = particles[j].mang[2]/jmax[2];
	      particles[j].je[0] = jmax[0];
	      particles[j].je[1] = jmax[1];
	      particles[j].je[2] = jmax[2];
	      /*if((particles[j].epsilon>9.0) && ((j<Header.Npart[0]) || (j>=Header.Npart[1]))) printf("PARTICULAS no circulares\n%lf %d %d %lf %lf %lf %lf\n",
												     particles[j].epsilon,j,particles[j].id,
												     particles[j].mang[2],
												     jmax,
												     sqrt(pow(particles[j].mang[0],2)+pow(particles[j].mang[1],2)+pow(particles[j].mang[2],2)),particles[j].Ebinding);*/
	    }
	}     
      if(jmax[2]!=-10.0)
	{
	  if(print==1)
	    {
	      fprintf(test1,"%e %e %e %e %e %e %e %e %d\n",e,jmax[2],Jmean,Jmean+2.5*sigma,Jmean-2.5*sigma,
		      particles[indice].pos[0],particles[indice].pos[1],particles[indice].pos[2],indice);
	    }
	}
      Emin = Emax;
    }
 
 
  /*
  if(print==1)
    {
      int counter;
      sprintf(namefiles,"epsilon_%d.dat",time);
      FILE *pf;
      pf = fopen(namefiles,"w");
      counter = 0;
      for(i=0; i<N_tot; i++)
	{
	  if(particles[i].type==2)
	    {
	      fprintf(pf,"%e %e %e %e\n",
		      particles[ind[counter]].Ebinding,
		      particles[ind[counter]].mang[2],
		      particles[ind[counter]].je[2],
		      particles[ind[counter]].epsilon);
	      counter = counter+1;
	    }
	  
	}
      
      fclose(pf);
      fclose(test);
      fclose(test1);
    }

*/
      
  free(r);
  free(p);
  free(E);
  free(jz);
  free(ind);

  return 0;
}

int Energybinding(int type,CM *centerm)
{

  size_t *IN,*p,ncore;
  int nparticles,i,j,k,*I;
  double *E,rcm[3],vcm[3],mass,v2;
  char outfile[1000];

  sprintf(outfile,"energy.%d",type);
  FILE *energy;
  energy = fopen(outfile,"w");
  
  
  nparticles = Header.npartTotal[type];
  ncore = 100;

  E = (double *)malloc((size_t)nparticles*sizeof(double));
  if(E == NULL){
    printf("Allocation error for E\n");
    exit(0);
  }

  IN = (size_t *)malloc((size_t)nparticles*sizeof(size_t));
  if(IN == NULL){
    printf("Allocation error for IN\n");
    exit(0);
  }

  p = (size_t *)malloc((size_t)nparticles*sizeof(size_t));
  if(p == NULL){
    printf("Allocation error for IN\n");
    exit(0);
  }

  I = (int *)malloc(ncore*sizeof(int));
  if(I == NULL){
    printf("Allocation error for I\n");
    exit(0);
  }

  
  k = 0;
  for(i=0; i<N_tot; i++)
    {
      if(particles[i].type==type)
	{
	  v2 = particles[i].vel[0]*particles[i].vel[0]
	    +particles[i].vel[1]*particles[i].vel[1]
	    +particles[i].vel[2]*particles[i].vel[2];

	  E[k] = 0.5*particles[i].mass*v2 + particles[i].mass*particles[i].pot;;
	  IN[k] = i;

	  fprintf(energy,"%e %e %e %e\n",
		  particles[i].pos[0],particles[i].pos[1],particles[i].pos[2],
		  E[k]); 
	  k = k+1;
	}
    }

  gsl_sort_index(p,E,1,nparticles);
  

  sprintf(outfile,"bound_particles.%d",type);
  FILE *binding;
  binding = fopen(outfile,"w");
  
  for(j=0; j<ncore; j++)
    {
      I[j] = IN[p[j]];
      fprintf(binding,"%d %e %e %e %e\n",
	      I[j],
	      particles[I[j]].pos[0],particles[I[j]].pos[1],particles[I[j]].pos[2],
	      E[p[j]]);
 
 
    }
  
  rcm[0]=rcm[1]=rcm[2]=0.0;
  vcm[0]=vcm[1]=vcm[2]=0.0;
  mass = 0.0;
  
  for(j=0; j<ncore; j++)
    {
      
      rcm[0] = rcm[0] + particles[I[j]].pos[0]*particles[I[j]].mass;
      rcm[1] = rcm[1] + particles[I[j]].pos[1]*particles[I[j]].mass;
      rcm[2] = rcm[2] + particles[I[j]].pos[2]*particles[I[j]].mass;
      
      vcm[0] = vcm[0] + particles[I[j]].vel[0]*particles[I[j]].mass;
      vcm[1] = vcm[1] + particles[I[j]].vel[1]*particles[I[j]].mass;
      vcm[2] = vcm[2] + particles[I[j]].vel[2]*particles[I[j]].mass;
      
      
      mass = mass + particles[I[j]].mass;	  
    }
  
  (*centerm).cm[0] = rcm[0]/mass;
  (*centerm).cm[1] = rcm[1]/mass;
  (*centerm).cm[2] = rcm[2]/mass;
  
  (*centerm).vcm[0] = vcm[0]/mass;
  (*centerm).vcm[1] = vcm[1]/mass;
  (*centerm).vcm[2] = vcm[2]/mass;
  

  printf("center type %d in ( %e , %e , %e )\n",type,
	 (*centerm).cm[0],
	 (*centerm).cm[1],
	 (*centerm).cm[2]);
  printf("   with velocity  ( %e , %e , %e )\n",
	 (*centerm).vcm[0],
	 (*centerm).vcm[1],
	 (*centerm).vcm[2]);
  
  
  fclose(energy);
  fclose(binding);
  free(E);
  free(p);
  free(I);
  free(IN);

  return 0;
}

int epsilonHistogram(char namefile[],int all,int type,int nbin)
{

  int i,j,N;
  double epsmax,epsmin,Max,Min,deltaeps;
  double histogram,ibin,fbin;
  char file[1000];
  
  //Finding min and max for epsilon
  epsmax = epsmin = 0.0;
  Max = -1e10;
  Min = 1e10;
  
  for(i=0; i<N_tot; i++) 
    {
      if((particles[i].type==2) || (particles[i].type==6))
	{
	  if((particles[i].r <= 5.0))
	    {    
	      if(particles[i].epsilon>Max)
		{
		  Max = particles[i].epsilon;
		  epsmax = particles[i].epsilon;
		}
	      else
		{
		  Max=Max;
		}
	      if(particles[i].epsilon<Min)
		{
		  Min = particles[i].epsilon;
		  epsmin = particles[i].epsilon;
		}
	      else
		{
		  Min=Min;
		}
	    }
	}
    }
  
  printf("epsilon min = %lf -> max = %lf )\n",epsmin,epsmax);
  
  deltaeps=(epsmax-epsmin)/nbin;
    
  ibin = epsmin;


  N = 0;
  for(i=0; i<N_tot; i++)
    {
      if(particles[i].type==2)
	{
	  N = N + 1;
	}
    }
  
  
  if(all==-1)
    {
      printf("\nHistogran epsilon with all partocles of disk\n");

      sprintf(file,"Histogram_eps_%s.dat",namefile);    
      FILE *his;
      his = fopen(file,"w");
      /*
      N = 0;
      for(i=0; i<N_tot; i++)
	{
	  if(particles[i].type==2)
	    {
	      N = N + 1;
	    }
	}
      */

      for(j=1; j<=nbin; j++)
	{
	  histogram = 0;
	  fbin = ibin + deltaeps;
	  
	  for(i=0; i<N_tot; i++)
	    {
	      if(particles[i].type==2)
		{
		  if((particles[i].epsilon>=ibin) && (particles[i].epsilon<fbin))
		    {
		      histogram = histogram + 1;
		    }
		}
	    }
	  
	  if((histogram>0))
	    {
	      fprintf(his,"%e %e %e\n",(ibin+fbin)/2.0,1.0*histogram/(1.0*N),sqrt(histogram));
	    }
	  ibin = fbin;
	}
      fclose(his);
     
    }
  else
    {
  
      printf("\nHistogran epsilon with particles of type %d\n",type);

      sprintf(file,"His_eps_%s_%d.dat",namefile,type);    
      FILE *his;
      his = fopen(file,"w");
      /*
      N = 0;
      for(i=0; i<N_tot; i++)
	{
	  if(particles[i].type==type)
	    {
	      N = N + 1;
	    }
	}
*/     
 
      for(j=1; j<=nbin; j++)
	{
	  histogram = 0;
	  fbin = ibin + deltaeps;
	  
	  for(i=0; i<N_tot; i++)
	    {
	      if(particles[i].type==type)
		{
		  if((particles[i].epsilon>=ibin) && (particles[i].epsilon<fbin))
		    {
		      histogram = histogram + 1;
		    }
		}
	    }
	  
	  if((histogram>0))
	    {
	      fprintf(his,"%e %e %e\n",(ibin+fbin)/2.0,1.0*histogram/(1.0*N),sqrt(histogram));
	    }
	  ibin = fbin;
	}
      fclose(his);
      
    }

  return 0;
}

//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
//GEOMETRY ROUTINES
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
/*
  ======================================================================
  Rotate the disk to the lagrangian plane
  ======================================================================
*/


/*
  ======================================================================
  convert cartesian to polar coodenates
  ======================================================================
*/


/*
  ======================================================================
  count the number of lines of a file
  ======================================================================
*/
int contadorParticulas(char *infile) 
{

  int NDAT,c;
  static FILE *pf;

  if((pf=fopen(infile,"r"))==NULL)
    printf("no puedo abrir archivo 1\n");

  NDAT=0;

  while((c=fgetc(pf))!= EOF)
    {
      if(c == '\n')
        {
          ++NDAT;
        }
    }

  fclose(pf);

  return NDAT;

}

int derivate(double *x, double *fx, double *dfdx, double *d2fdx2, int npuntos)
{
  double dx;
  int i;

  for(i=2; i<npuntos-2; i++)
    {
      //intervalos
      dx = x[i+1]-x[i];
      //Primera derivada
      dfdx[i-2] = (fx[i-2] -8.0*fx[i-1] +8.0*fx[i+1] -fx[i+2])/(12.0*dx);
      //Segunda derivada
      d2fdx2[i-2] = (-fx[i+2] +16.0*fx[i+1] -30.0*fx[i] +16.0*fx[i-1] -fx[i-2])/(12.0*dx*dx);
    }
  printf("derivada normal\n");

  //liberar dfdx y d2fdx2 en el codigo donde se utilice esta rutina

  return 0;
}

int derivateLog(double *x, double *fx, double *dfdx, double *d2fdx2, int npuntos)
{
  int i;
  double dx;
 
  for(i=2; i<npuntos-2; i++)
    {
      //intervalos
      dx = log10(x[i+1]) - log10(x[i]);
      //Primera derivada
      dfdx[i-2] = (fx[i-2] -8.0*fx[i-1] +8.0*fx[i+1] -fx[i+2])/(12.0*dx);
      //Segunda derivada
      d2fdx2[i-2] = (-fx[i+2] +16.0*fx[i+1] -30.0*fx[i] +16.0*fx[i-1] -fx[i-2])/(12.0*dx*dx);
    }
  printf("derivada logaritmica\n");
  return 0;
}

int gslSmooth(double *xi, double *fxi, int npuntos)
{
  size_t k;
  k = 4;
  size_t n = npuntos;
  size_t ncoeffs = 12;
  size_t nbreak = ncoeffs +2 -k;
  size_t i, j;
  gsl_bspline_workspace *bw;
  gsl_vector *B;
  //  double dy;
  //gsl_rng *r;
  gsl_vector *c, *w;
  gsl_vector *x, *y;
  gsl_matrix *X, *cov;
  gsl_multifit_linear_workspace *mw;
  double chisq, Rsq, dof, tss;
  
  //gsl_rng_env_setup();
  //r = gsl_rng_alloc(gsl_rng_default);
  printf("%zd %zd %zd\n",n,ncoeffs,nbreak);
  /* allocate a cubic bspline workspace (k = 4) */
  bw = gsl_bspline_alloc(4, nbreak);
  B = gsl_vector_alloc(ncoeffs);
  x = gsl_vector_alloc(n);
  y = gsl_vector_alloc(n);
  X = gsl_matrix_alloc(n, ncoeffs);
  c = gsl_vector_alloc(ncoeffs);
  w = gsl_vector_alloc(n);
  cov = gsl_matrix_alloc(ncoeffs, ncoeffs);
  mw = gsl_multifit_linear_alloc(n, ncoeffs);
  
  //printf("#m=0,S=0\n");
  /* this is the data to be fitted */
  for (i=0; i<n; ++i)
    {
      double sigma;
      //double xi = (15.0 / (N - 1)) * i;
      //double yi = cos(xi) * exp(-0.1 * xi);
      
      sigma = 0.1 * fxi[i];
      //dy = gsl_ran_gaussian(r, sigma);
      //yi += dy;
      
      gsl_vector_set(x, i, xi[i]);
      gsl_vector_set(y, i, fxi[i]);
      gsl_vector_set(w, i, 1.0 / (sigma * sigma));
      
      //    printf("%lf %lf\n", xi[i], fxi[i]);
    }
  
  /* use uniform breakpoints on [0, 15] */
  gsl_bspline_knots_uniform(0.0, 10.0, bw);
  
  /* construct the fit matrix X */
  for (i=0; i<n; ++i)
    {
      double xi = gsl_vector_get(x, i);
      
      /* compute B_j(xi) for all j */
      gsl_bspline_eval(xi, B, bw);
      
      /* fill in row i of X */
      for (j = 0; j < ncoeffs; ++j)
	{
	  double Bj = gsl_vector_get(B, j);
	  gsl_matrix_set(X, i, j, Bj);
	}
    }
  
  /* do the fit */
  gsl_multifit_wlinear(X, w, y, c, cov, &chisq, mw);
  
  dof = n - ncoeffs;
  tss = gsl_stats_wtss(w->data, 1, y->data, 1, y->size);
  Rsq = 1.0 - chisq / tss;
  
  fprintf(stderr, "chisq/dof = %e, Rsq = %f\n", 
	  chisq / dof, Rsq);
  
  /* output the smoothed curve */
  {
    double yi, yerr;
    //double xi, yi, yerr;
    
    //printf("#m=1,S=0\n");
    FILE *salida;
    salida = fopen("pot_suavisado.dat","w");
    //for (xi = 0.0; xi < 10.0; xi += 0.1)
    for(i=0; i<npuntos; i++)
      {
	gsl_bspline_eval(xi[i], B, bw);
	//gsl_bspline_eval(xi, B, bw);
	gsl_multifit_linear_est(B, c, cov, &yi, &yerr);
	fprintf(salida,"%lf %lf\n", xi[i], yi);
      }
    fclose(salida);
  }
  
  //gsl_rng_free(r);
  gsl_bspline_free(bw);
  gsl_vector_free(B);
  gsl_vector_free(x);
  gsl_vector_free(y);
  gsl_matrix_free(X);
  gsl_vector_free(c);
  gsl_vector_free(w);
  gsl_matrix_free(cov);
  gsl_multifit_linear_free(mw);
  
  return 0;
}
/*
int derivateLogjuank(double *x, double *fx, int npuntos)
{
  int i;
  double dx;

  //intervalos
  dx = log10(x[2]) - log10(x[1]);
  
  for(i=0; i<npuntos; i++)
    {
      
      // Primera derivada
      if(i < (npuntos-2))
	dfdx[i] = (fx[i+1] -fx[i])/dx;
      if(i == (npuntos-1))
	dfdx[i] = (fx[i] -fx[i])/dx;
      
      //Segunda derivada
      if( (i >= 1) && (i < (npuntos-2)) )
	d2fdx2[i] = (fx[i+1] - 2.0*fx[i] + fx[i-1])/(dx*dx);
      if( (i == 0) )
	d2fdx2[i] = (fx[i+1] - 2.0*fx[i] + fx[i])/(dx*dx);
      if( (i == (npuntos-1)) )
	d2fdx2[i] = (fx[i] - 2.0*fx[i] + fx[i])/(dx*dx);
    }
  
  printf("derivada logaritmica\n");
  return 0;
}
*/
int totalmass(char *infile)
{
  int NTOTAL=0,npart[6],n,i,j;
  double mass[6],masstotal=0.0;
  particulas *diskpart;
  char infiles[1000];

  readGadgetBinary("option_makefile.dat","small_galaxy_000");

  for(i=0; i<6; i++)
    {
      NTOTAL = NTOTAL+Header.npartTotal[i];
    }

  diskpart = (particulas *)malloc((size_t)NTOTAL*sizeof(particulas));
  if(diskpart == NULL){
    printf("Allocation error routines:81\n");
    exit(0);
  }

  n=0;
  for(j=0; j<6; j++)
    {
      npart[j] = 0;
      mass[j] = 0.0;
      if( (Header.npartTotal[j])!=0 )
	{
	  FILE *diskfile;
	  sprintf(infiles,"%s.%d",infile,j);
	  diskfile = fopen(infiles,"r");
	  npart[j] = contadorParticulas(infiles);
	  
	  for(i=n; i<n+npart[j] ; i++)
	    {
	      fscanf(diskfile,"%d %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e\n",
		     &diskpart[i].id,
		     &diskpart[i].pos[0],&diskpart[i].pos[1],&diskpart[i].pos[2],
		     &diskpart[i].vel[0],&diskpart[i].vel[1],&diskpart[i].vel[2],
		     &diskpart[i].mass,
		     &diskpart[i].pot,
		     &diskpart[i].posc[0],&diskpart[i].posc[1],&diskpart[i].posc[2],
		     &diskpart[i].velc[0],&diskpart[i].velc[1],&diskpart[i].velc[2],
		     &diskpart[i].mang[0],&diskpart[i].mang[1],&diskpart[i].mang[2],
		     &diskpart[i].mr,
		     &diskpart[i].r,
		     &diskpart[i].vc,
		     &diskpart[i].epsilon);
	      mass[j] = mass[j]+diskpart[i].mass;
	    }
	  fclose(diskfile);
	}
      masstotal = masstotal+mass[j];
      n = n + npart[j];
      
      printf("tipo %d with %d particles of mass %e = %e\n",j,npart[j],Header.mass[j],mass[j]);
    
    }

  return 0;
}


int rotation_curve(int type,int Npart, char outfile[])
{
  int i,counter;

  double *r,*mr,vc,m;
  size_t *I,*P;
  
  char namefiles[1000];

  FILE *Mr;

  sprintf(namefiles,"rotation_curve_%s.%d",outfile,type);
  Mr = fopen(namefiles,"w");
   
  printf("%d\t",type);
  printf("%d\n",Npart);
  
  r = (double *) malloc(Npart*sizeof(double));
  if(r == NULL) 
    {
      printf("No memory for r\n"); 
      exit(0);
    }
  
  mr = (double *) malloc(Npart*sizeof(double));
  if(mr == NULL) 
    {
      printf("No memory for mr\n"); 
      exit(0);
    }
  
  I = (size_t *)malloc(Npart*sizeof(size_t));
  if(I == NULL) 
    {
      printf("No memory for I\n");
      exit(0);
    }
  
  P = (size_t *)malloc(Npart*sizeof(size_t));
  if(P == NULL) 
    {
      printf("No memory for P\n");
      exit(0);
    }
  
  counter = 0;
  for(i=0; i<N_tot; i++)
    {
      if(particles[i].type==type)
	{
	  r[counter] = particles[i].r;
	  I[counter] = i;
	  counter++;
	}
    }
  
  gsl_sort_index(P,r,1,Npart);
   
  particles[I[P[0]]].mr = 0.0;
  m = particles[I[P[0]]].mass;
  counter = 0;
  for(i=0; i<N_tot; i++)
    {
      if(particles[i].type==type)
	{
	  particles[I[P[counter]]].mr = m;
	  m = m + particles[I[P[counter]]].mass;
	  counter++;
	}
    }
  
  counter = 0;
  for(i=0; i<N_tot; i++)
    {
      if(particles[i].type==type)
	{
	  mr[counter] = particles[I[P[counter]]].mr;
	  vc = GGAD*mr[counter]/r[P[counter]];
	  fprintf(Mr,"%e %e %e\n",r[P[counter]],mr[counter],vc);
	  counter++;
	}
    }
  
  free(r);
  free(mr);
  free(I);
  free(P);
  fclose(Mr);
   
  return 0;
}

int total_rotation_curve(char outfile[])
{
  int i;

  double *r,*mr,vc,m;
  size_t *I,*P;
  
  char namefiles[1000];

  FILE *Mr;

  sprintf(namefiles,"rotation_curve_%s_total.dat",outfile);
  Mr = fopen(namefiles,"w");
     
  r = (double *) malloc(N_tot*sizeof(double));
  if(r == NULL) 
    {
      printf("No memory for r\n"); 
      exit(0);
    }
  
  mr = (double *) malloc(N_tot*sizeof(double));
  if(mr == NULL) 
    {
      printf("No memory for mr\n"); 
      exit(0);
    }
  
  I = (size_t *)malloc(N_tot*sizeof(size_t));
  if(I == NULL) 
    {
      printf("No memory for I\n");
      exit(0);
    }
  
  P = (size_t *)malloc(N_tot*sizeof(size_t));
  if(P == NULL) 
    {
      printf("No memory for P\n");
      exit(0);
    }
  
  for(i=0; i<N_tot; i++)
    {
      r[i] = particles[i].r;
      I[i] = i;
    }
  
  gsl_sort_index(P,r,1,N_tot);
   
  particles[I[P[0]]].mr = 0.0;
  m = particles[I[P[0]]].mass;

  for(i=0; i<N_tot; i++)
    {
      particles[I[P[i]]].mr = m;
      m = m + particles[I[P[i]]].mass;
    }
  
  for(i=0; i<N_tot; i++)
    {
      mr[i] = particles[I[P[i]]].mr;
      vc = GGAD*mr[i]/r[P[i]];
      fprintf(Mr,"%e %e %e\n",r[P[i]],mr[i],vc);
    }
  
  free(r);
  free(mr);
  free(I);
  free(P);
  fclose(Mr);
   
  return 0;
}

//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
//MORPHOLOGICAL ROUTINES
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

/*
  Usage:
      readGadgetParam("MW_CLUES_10909_1_MMP_WB.param","SofteningDisk",'f',&T);
      readGadgetParam("MW_CLUES_10909_1_MMP_WB.param","ICFormat",'d',&d);
      readGadgetParam("MW_CLUES_10909_1_MMP_WB.param","InitCondFile",'s',cadena);
 */
/*
void readGadgetParam(char *paramfile,char *variable,char type,void *result)
{
  char tmpfile[]="/tmp/readgad";
  char cmd[1000];
  char value[100];
  //float tmp;

  sprintf(cmd,"grep '^%s ' %s | awk '{print $2}' > %s",variable,paramfile,tmpfile);
  system(cmd);
  FILE *fp;
  fp=fopen(tmpfile,"r");
  fscanf(fp,"%s",value);
  fclose(fp);
  fprintf(stdout,"Value read from param. file : %s\n",value);
  switch(type){
  case 'f':{
    float *valread=(float*)result;
    *valread=atof(value);
    break;
  }
  case 'd':{
    int *valread=(int*)result;
    *valread=atoi(value);
    break;
  }
  case 's':{
    char *valread=(char*)result;
    strcpy(valread,value);
    break;
  }
  }
}
*/
void logBins(float xmin,float xmax,int nbin,float *dlogx)
{
  *dlogx=(log10(xmax)-log10(xmin))/nbin;
}

float gravSoft(float x,float h)
{
  
  float kern=0.0,u;

  u=x/h;

  if( (u >= 0) && (u < 0.5) )
    kern = 
      (16.0/3.0)*u*u - 
      (48.0/5.0)*pow(u,4) + 
      (32.0/5.0)*pow(u,5) - (14.0/5.0); 
  else if((u >= 0.5) && (u < 1.0))
    kern= 
      (1.0/(15.0*u)) + 
      (32.0/3.0)*u*u - 
      16.0*pow(u,3) + 
      (48.0/5.0)*pow(u,4) - 
      (32.0/15.0)*pow(u,5) - 
      (16.0/5.0);
  else if(u >= 1.0)
    kern=-1.0/u;
  
  return (-kern);
}

float distanceVectors(vector v1,vector v2)
{
  float dx,dy,dz;
  dx=v1.comp[0]-v2.comp[0];
  dy=v1.comp[1]-v2.comp[1];
  dz=v1.comp[2]-v2.comp[2];
  //fprintf(stdout,"%e,%e,%e\n",dx,dy,dz);
  return sqrt(dx*dx+dy*dy+dz*dz);
}

int dispersionVelocity(int nbines)
{
  double *R,sigmaVR,meanVR,delta;
  double Rmin,Rmax,Rini,Rfin;
  int i,j,npart;
  size_t *indexR;
  FILE *disvel;

  nbines=10;

  R=(double*)malloc(N_part[2]*sizeof(double));
  indexR=(size_t*)malloc(N_part[2]*sizeof(size_t));
  
  FILE *vel;
  vel=fopen("velocidades_disco","w");
  for(i=0;i<N_part[2];i++)
    {
      R[i]=Pos[i+N_part[1]].comp[3];
      indexR[i]=Id[i+N_part[1]];
      //p int velocities of disk
         
      fprintf(vel,"%lf %lf %lf %lf %lf %lf\n",
	      Vel[i+N_part[1]].comp[0],
	      Vel[i+N_part[1]].comp[1],
	      Vel[i+N_part[1]].comp[2],
	      Vel[i+N_part[1]].comp[3],
	      Vel[i+N_part[1]].comp[4],
	      Vel[i+N_part[1]].comp[5]);
    }

  fclose(vel);
  
  gsl_sort_index(indexR,R,1,N_part[2]);
  
  //Rmin=log10(R[indexR[0]]);
  // Rmax=log10(R[indexR[N_part[2]-1]]);
  //delta=(Rmax-Rmin)/nbines;

  Rmin=R[indexR[0]];
  Rmax=R[indexR[N_part[2]-1]];
  delta=(Rmax-Rmin)/nbines;

  printf("\n numero de bines %d\n",nbines);

  disvel=fopen("dispersion_velocities_disk.dat","w");
  Rini=Rmin;
  for(j=0;j<nbines;j++)
    {
      Rfin=Rini+delta;
      meanVR=0.0;
      sigmaVR=0.0;
      npart=0.0;
      for(i=0;i<N_part[2];i++)
	{
	  //if((R[indexR[i]]>=pow(10,Rini)) && (R[indexR[i]]<=pow(10,Rfin)))
	  if((R[indexR[i]]>=Rini) && (R[indexR[i]]<=Rfin))
	    {
	      meanVR=meanVR+Vel[indexR[i]].comp[3];
	      npart=npart+1;
	    }    	  
	}
      meanVR=meanVR/npart;
      npart=0;
      for(i=0;i<N_part[2];i++)
	{
	  //if((R[indexR[i]]>=pow(10,Rmin)) && (R[indexR[i]]<=pow(10,Rfin)))
	  if((R[indexR[i]]>=Rmin) && (R[indexR[i]]<=Rfin))
	    {
	      sigmaVR=sigmaVR+(meanVR-Vel[indexR[i]].comp[3]);
	      npart=npart+1;
	    }
	}
     sigmaVR = sigmaVR/npart;
      if(npart>0)
	{
	  // fprintf(disvel,"%lf %lf %lf\n",pow(10,Rini),sigmaVR,meanVR);
	  fprintf(disvel,"%lf %lf %lf\n",Rini,log10(sigmaVR),meanVR);
	}
      Rini=Rfin;
    }
  fclose(disvel);
  return 0;
}

int diskPotential(float rmin,float rmax,int nbines)
{
  FILE *fp=fopen("radial-potential.dat","w");
  int i,m,k,n,nphisamples,nzsamples;
  float dlogr,logrb,logre;
  float rbin,phi,zmin,zmax,h;
  float phim,dist,phiavg;
  vector rsample;

  int npa,npb;
  float npe;

  npa=0;
  npb=3;
  npe=1./3;

  fprintf(stdout,"Computing the potential at the disk plane in %d bins...\n",nbines);

  if(!NPotential){
    Potential=(float*)malloc(nbines*sizeof(float));
    RPotential=(float*)malloc(nbines*sizeof(float));
    NPotential=nbines;
  }

  dlogr=(log10(rmax)-log10(rmin))/nbines;

  //fprintf(stdout,"rmin,rmax,dlogr=%e,%e,%e\n",rmin,rmax,dlogr);

  //LOOP ON THE BINS
  logrb=log10(rmin);
  for(i=0;i<nbines;i++){
    logre=logrb+dlogr;
    rbin=pow(10,(logre+logrb)/2);

    //LOOP ON SAMPLE POINTS
    phiavg=0;
    nphisamples=(int)(npa*pow((rbin-rmin)/rmin,npe)+npb);
    nzsamples=1;
    zmin=0.0;zmax=0.0;

    //fprintf(stdout,"Bin %d (n.phi.samples=%d): [%e,%e] (R = %e)\n",i,nphisamples,logrb,logre,rbin);
    for(m=0;m<nphisamples;m++){
      phi=2*M_PI/nphisamples*m;
      rsample.comp[0]=rbin*cos(phi);
      rsample.comp[1]=rbin*sin(phi);

      for(k=0;k<nzsamples;k++){
	rsample.comp[2]=k*(zmax-zmin)/nzsamples+zmin;
	phim=0;
	
	//fprintf(stdout,"\tSample point (%d,%d): phi = %e, z = %e\n",m,k,phi,rsample.comp[2]);
	for(n=0;n<N_tot;n++){
	  dist=distanceVectors(Pos[n],rsample);
	  /*
	  fprintf(stdout,"\t\tParticle %d: r = (%e,%e,%e)\n",
		  n,Pos[n].comp[0],Pos[n].comp[1],Pos[n].comp[2]);
	  fprintf(stdout,"\t\tDistance between %d and sample point: %e\n",n,dist);
	  */
	  h=Softening[n];
	  phim-=GGAD*Mass[n]*gravSoft(dist,h)/h;
	}
	//fprintf(stdout,"\tTotal potential in sample point: phim = %e\n",phim);
	phiavg+=phim;
      }
    }
    phiavg/=(nphisamples*nzsamples);
    //fprintf(stdout,"\tAverage potential in bin %d (n=%d)= %e\n",i,nphisamples*nzsamples,phiavg);
    logrb=logre;
    RPotential[i]=rbin;
    Potential[i]=phiavg;
    fprintf(fp,"%e %e\n",rbin,phiavg);
    //exit(0);
  }
  fclose(fp);
  return 0;
}
int printType(char *inputfile,int type)
{
  int i;
  char outfiles[1000];

  sprintf(outfiles,"%s",inputfile);
  printf("Printing particles type %d in %s\n",type,outfiles);
  FILE *f;
  f = fopen(outfiles,"w");

  for(i=0; i<N_tot; i++)
    {
      if(particles[i].type==type)
	{
	  fprintf(f,"%d %e %e %e %e %e %e %e %e\n",
		  particles[i].id,
		  particles[i].pos[0],particles[i].pos[1],particles[i].pos[2],
		  particles[i].vel[0],particles[i].vel[1],particles[i].vel[2],
		  particles[i].mass,
		  particles[i].pot);
	}
    }
  fclose(f);

 return 0;
}
//
int printAscii(char *inputfile,int all,int type)
{
  int i,j,ini,fin;
  char outfiles[100];
  
  if(type!=-1)
    {
      sprintf(outfiles,"%s.%d",inputfile,type);
      FILE *f;
      f = fopen(outfiles,"w");
      
      for(i=0; i<N_tot; i++)
	{
	  if(particles[i].type==type)
	    {
	      fprintf(f,"%d %e %e %e %e %e %e %e %e\n",
		      particles[i].id,
		      particles[i].pos[0],particles[i].pos[1],particles[i].pos[2],
		      particles[i].vel[0],particles[i].vel[1],particles[i].vel[2],
		      particles[i].mass,
		      particles[i].pot); 
	    }
	  
	}
      fclose(f);
    }
  else
    {
      ini = 0;
      fin = 0;
      if(all==0)
	{
	  for(i=0; i<6; i++)
	    {
	      fin = fin + Header.npartTotal[i];
	      if(Header.npartTotal[i]!=0)
		{
		  sprintf(outfiles,"%s.%d",inputfile,i);
		  printf("Printing %d particles in %s ",Header.npartTotal[i],outfiles);
		  FILE *f;
		  f = fopen(outfiles,"w");
		  printf("with id's between %d - %d\n",ini,fin);
		  for(j=ini; j<fin; j++)
		    {
		      fprintf(f,"%d %e %e %e %e %e %e %e %e\n",
			      particles[j].id,
			      particles[j].pos[0],particles[j].pos[1],particles[j].pos[2],
			      particles[j].vel[0],particles[j].vel[1],particles[j].vel[2],
			      particles[j].mass,
			      particles[j].pot);
		    }
		  fclose(f);
		}
	      ini = fin;
	    }
	}
      if(all==1)
	{
	  for(i=0; i<=6; i++)
	    {
	      fin = fin + Header.npartTotal[i];
	      
	      //if(Header.npartTotal[i]!=0)
	      //{
	      sprintf(outfiles,"%s_all.%d",inputfile,i);
	      printf("Printing %d particles in %s with all properties",Header.npartTotal[i],outfiles);
	      FILE *f;
	      f = fopen(outfiles,"w");
	      printf("with id's between %d - %d\n",ini,fin);
	      for(j=0; j<N_tot; j++)
		{
		  if(particles[j].type==i)
		    {
		      fprintf(f,"%d %d %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e\n",
			      particles[j].id,
			      particles[j].type,
			      particles[j].pos[0],particles[j].pos[1],particles[j].pos[2],
			      particles[j].vel[0],particles[j].vel[1],particles[j].vel[2],
			      particles[j].mass,
			      particles[j].pot,
			      particles[j].posc[0],particles[j].posc[1],particles[j].posc[2],
			      particles[j].velc[0],particles[j].velc[1],particles[j].velc[2],
			      particles[j].mang[0],particles[j].mang[1],particles[j].mang[2],
			      particles[j].mr,
			      particles[j].r,
			      particles[j].vc,
			      particles[j].epsilon,
			      particles[j].Ebinding,
			      particles[j].je[0],particles[j].je[1],particles[j].je[2]);
		    }
		}
	      fclose(f);
	      //}
	      ini = fin;
	    }
	}
    }
  return 0;
}

int inclination()
{

  int i,n,counter,Izmin;
  size_t *p;
  double Md,Mb,h,zo,a;
  double A,B,C,x,y,z,R,zmin,sigma,vzmax;
  double alpha1,alpha2,nu,vz,Lz,L;
  double *Z;

  FILE *pf;
  pf = fopen("inclinations.dat","w");

  n = 0; 
  for(i=0; i<N_tot; i++)
    {
      if(particles[i].type==2)
	{
	  n = n+1;
	}
    }

  p = (size_t *)malloc(n*sizeof(size_t));
  if(p == NULL) 
    {
      printf("No memory for p\n");
      exit(0);
    }
  Z = (double *)malloc(n*sizeof(double));
  if(Z == NULL) 
    {
      printf("No memory for Z\n");
      exit(0);
    }
 
  counter = 0;
  Md = 0.0;
  for(i=0; i<N_tot; i++)
    {
      if(particles[i].type==2)
	{
	  Z[counter] = fabs(particles[i].pos[2]);
	  Md = Md + particles[i].mass; 
	  counter++;
	}
    }
  

  gsl_sort_index(p,Z,1,n);


  Mb = 0.2*Md;
  h = 1.84438;
  zo = 0.2*h;
  a = 0.2*h; 

  zmin = zo;
  vzmax = -1.0; 
  for(i=0; i<N_tot; i++)
    {
      if(particles[i].type==2)
	{
	  if((particles[i].posc[0]>(2.0*h-0.01)) && (particles[i].posc[0]<(2.0*h+0.01)))
	    {
	      if(fabs(particles[i].vel[2])>vzmax)
		{
		  if(fabs(particles[i].pos[2])<0.1*zo)
		    {
		      zmin = fabs(particles[i].pos[2]);
		      vzmax = particles[i].vel[2];
		      Izmin = i;
		    }
		  else
		    {
		      zmin = zmin;
		      vzmax = vzmax;
		    }
		}
	     
	    }
	}
    }

 

  for(i=0; i<N_tot; i++)
    {
      if(particles[i].type==2)
	{
	  // if((particles[i].posc[0]>(2.0*h-0.01)) && (particles[i].posc[0]<(2.0*h+0.01)))
	  //{
	      vz = particles[i].vel[2];
	      sigma = vzmax*vzmax-vz*vz;

	      R = particles[i].posc[0];
	      x = particles[i].pos[0];
	      y = particles[i].pos[1];
	      z = particles[i].pos[2];

	      A = 1.0/(2.0*h*h*zo*z);
	      B = 0.2/sqrt(R*R+z*z);
	      C = 1.0/pow((sqrt(R*R+z*z)+a),2);

	      nu = A*exp(-R/h)+B*C*GGAD*Md;

	      Lz = particles[i].mang[2];
	      L = sqrt(particles[i].mang[0]*particles[i].mang[0]+
		       particles[i].mang[1]*particles[i].mang[1]+
		       particles[i].mang[2]*particles[i].mang[2]);
	     
	      alpha1 = (M_PI/2.0)-acos(fabs(Lz)/L);

	      alpha2 = atan(sqrt(sigma)/(R*nu));

	      if(sigma>=0.0)
		{
		  fprintf(pf,"%e %e %e %e %e %e %e %e %e %e\n",x,y,z,R,vz,fabs(vzmax)-fabs(vz),nu,sqrt(sigma),alpha1,alpha2);
		}
	      //}
	}
    }


  /*
  counter = 0;
  for(i=0; i<N_tot; i++)
    {
      if(particles[i].type==2)
	{
	  R = particles[p[counter]].posc[0];
	  z = Z[p[counter]];

	  A = 1.0/(2.0*h*h*zo*z);
	  B = 0.2/sqrt(R*R+z*z);
	  C = 1.0/pow((sqrt(R*R+z*z)+a),2);


	  nu = A*exp(-R/h)+B*C*GGAD*Md;
 
	  vz = particles[p[counter]].vel[2];
	  E = fabs(0.5*vz*vz + particles[p[counter]].pot);
	  
	  //E = fabs(particles[p[counter]].Ebinding);
	 
	  f1 = 1.0/(R*nu);
	  f2 = 2.0*E*nu-vz*vz;
	  //	  alpha = atan(f1/sqrt(f2));

	  Lz = particles[i].mang[02];
	  L = sqrt(particles[i].mang[0]*particles[i].mang[0]+
		   particles[i].mang[1]*particles[i].mang[1]+
		   particles[i].mang[2]*particles[i].mang[2]);
	  alpha = (M_PI/2.0)-acos(Lz/L);

	  fprintf(pf,"%e %e %e %e %e %e %e %e\n",Z[p[counter]],R,vz,nu,E,E*nu,f2,alpha);
	  counter++;
	}
    }
*/

  free(p);
  free(Z);
  fclose(pf);

  return 0;
}
